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This  paper  presents  a  degradation-based  procedure  for  estimating  the  full  and  residual  lifetime  distribution  of  a  single-unit  system 
subject  to  Markovian  deterioration.  The  hybrid  approach  unites  real  degradation  measures  with  analytical,  stochastic  failure  models 
to  numerically  compute  the  distributions  and  their  moments.  Two  distinct  models  are  shown  to  perform  well  when  compared  with 
simulated  data.  Moreover,  results  obtained  from  the  second  model  are  compared  with  empirically  generated  lifetimes  of  2024- 
T3  aluminum  alloy  specimens.  The  numerical  experiments  indicate  that  the  proposed  techniques  are  useful  for  remaining  lifetime 
prognosis  in  both  cases. 


1.  Introduction 

In  this  paper,  we  consider  the  problem  of  using  sensor  data 
to  estimate  full  and  residual  lifetime  distributions  for  a 
single-unit  system  subject  to  a  stochastically  evolving  envi¬ 
ronment.  We  are  particularly  motivated  by  applications  in 
autonomous  prognostic  systems  whose  primary  goals  are 
to  assess  the  current  and  future  health  of  single-  or  multi¬ 
unit  systems  to  avoid  catastrophic  failures  and  to  elimi¬ 
nate  superfluous  preventive  maintenance  activities.  It  is  well 
known  that  the  difficulty  in  estimating  lifetime  probabilities 
is  the  scarcity  or  absence  of  actual  failure  time  observations 
(failure-based  reliability).  This  stems  from  the  fact  that  it  is 
often  infeasible  to  run  components  to  failure,  whereas  ac¬ 
celerated  life  tests  may  not  be  representative  of  the  true  op¬ 
erating  environment  of  components;  thus,  yielding  possibly 
unreliable  results.  Hence,  novel  techniques  that  exploit  read¬ 
ily  available  sensor  data  for  lifetime  estimation  are  needed. 
Examples  of  such  data  may  be  the  current  state  of  the  am¬ 
bient  environment  of  the  single-unit  system  or  measures  of 
degradation  suffered  by  the  single-unit  system  up  to  some 
point  in  time. 

Relative  to  failure-based  reliability,  degradation-based 
reliability  has  received  a  modest  amount  of  attention  in 
the  open  literature.  The  use  of  degradation  measures  to  as¬ 
sess  component  lifetimes  was  addressed  in  the  early  work  of 
Gertsbakh  and  Kordonsky  (1969)  who  used  a  sample  path 
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approach  to  assess  the  reliability  for  a  simple  linear  degra¬ 
dation  path  with  random  slope  and  intercept.  Later,  Lu  and 
Meeker  (1993)  estimated  the  failure  time  distribution  of  a 
unit  that  fails  once  the  cumulative  damage  exceeds  a  given 
threshold  value.  Their  parametric  statistical  model  allows 
for  a  bootstrap  method  to  generate  confidence  intervals 
around  the  failure  time  distribution,  but  does  not  necessar¬ 
ily  allow  for  a  closed-form  expression  of  the  distribution. 
More  recently.  Meeker  and  Escobar  (1998)  provide  a  use¬ 
ful  summary  of  degradation  models,  emphasizing  the  use 
of  linear  models  with  assumed  log-normal  rates  of  degra¬ 
dation.  In  such  a  case,  the  full  lifetime  distribution  can  be 
computed  analytically.  Other  recent  models  encountered  in 
the  literature  deal  with  the  degradation  of  materials  such  as 
those  due  to  Gillen  and  Celina  (2001).  Meeker  et  al.  (2002) 
discussed  general  approaches  to  estimating  lifetime  distri¬ 
butions  in  accelerated  life  tests  for  highly  variable  environ¬ 
ments.  The  models  presented  therein  focus  on  the  specifi¬ 
cation  of  the  degradation  path  as  it  depends  explicitly  on 
the  operating  environment.  In  that  paper,  the  authors  note 
the  need  for  a  formal  stochastic  process  model  of  the  am¬ 
bient  environment  and  for  expedient  numerical  techniques 
for  the  evaluation  of  these  measures.  This  work  provides  a 
first  step  in  addressing  both  of  these  issues. 

Analytical  lifetime  distribution  models  for  single-unit 
systems  have  been  studied  extensively  and  exist  primarily  in 
the  form  of  stochastic  shock  and  wear  models.  Bogdanoff 
and  Kozin  (1985)  provide  a  good  summary  of  probabilistic 
models  of  cumulative  damage  focusing  on  discrete-time  ver¬ 
sions  of  shock  models.  Such  models  assume  that  the  unit 
sustains  a  random  amount  of  damage  each  time  a  shock 
occurs  at  either  fixed  or  random  time  intervals.  Wear  pro¬ 
cesses  differ  from  shock  models  in  that  they  assume  that 
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damage  accumulates  continuously  over  time.  Some  prop¬ 
erties  of  the  failure  time  distribution  have  been  examined 
for  such  processes  when  the  degradation  path  is  assumed 
to  be  a  Levy  process  (Esary  et  al.,  1973;  Abdel-Elameed, 
1984).  Singpurwalla  (1995)  gave  an  excellent  summary  of 
a  variety  of  stochastic  failure  time  models  for  systems  in  a 
random  environment  and  particularly  noted  the  difficulty 
of  implementing  many  of  these  in  a  practical  setting.  Re¬ 
cently,  Kharoufeh  (2003)  provided  explicit  transform  re¬ 
sults  for  the  full  lifetime  distribution  and  moments  for  a 
single-unit  system  in  a  time-homogeneous  Markovian  en¬ 
vironment.  The  main  results  of  that  work  provide  the  basis 
for  our  procedures  here. 

In  the  spirit  of  other  degradation-based  approaches,  we 
assume  that  a  single-unit  system  resides  in  a  random  oper¬ 
ating  environment  so  that  its  rate  of  degradation  depends 
explicitly  on  the  influence  of  external  factors.  The  evolution 
of  the  random  environment  is  characterized  as  a  stationary 
Continuous-Time  Markov  Chain  (CTMC).  The  key  inno¬ 
vation  of  our  approach  is  the  union  of  analytical  models 
(specifically  Markovian  wear  processes)  with  the  use  of  real 
sensor  data  to  estimate  full  and  residual  lifetime  distribu¬ 
tions.  We  present  two  models.  The  first  one  assumes  that 
the  sensor  data  provide  information  on  the  state  of  the  ran¬ 
dom  environment.  The  second  one  provides  information  on 
the  cumulative  degradation  up  to  some  fixed  point  in  time. 
Our  aim  is  to  provide  a  novel,  hybrid  approach  that  will  ul¬ 
timately  lead  to  a  feasible,  and  implementable,  prognostic 
system. 

The  remainder  of  the  paper  is  organized  as  follows.  The 
next  section  reviews  the  mathematical  model  describing 
unit  lifetime  distributions.  In  Section  3  we  demonstrate  the 
means  by  which  real  sensor  data  may  be  used  to  estimate 
the  model  parameters  when  the  operating  environment  is 
observable.  Section  4  presents  a  model  for  the  case  when 
only  the  cumulative  degradation  level  is  observable.  Several 
illustrative  examples  are  provided  in  Section  5,  and  our  con¬ 
cluding  remarks  and  future  research  directions  are  given  in 
Section  6. 


2.  Mathematical  model  description 

This  section  provides  the  general  mathematical  framework 
upon  which  we  build  our  estimation  technique  using  real 
degradation  data.  It  briefly  summarizes  the  main  results 
of  Kharoufeh  (2003)  for  a  single-unit  system  subject  to  a 
random  environment  modelled  as  a  CTMC  and  extends  the 
results  to  the  residual  lifetime  of  the  unit. 

The  rate  of  degradation  of  the  system  at  time  t  >  0  is 
governed  by  a  random  environment  that  is  modelled  as 
a  finite  state  CTMC,  {Z(t)  :  t  >  0}.  Let  Z(t)  denote  the 
state  of  the  random  environment  process  at  time  t  e  M+, 
where  M ,  denotes  the  non-negative  real  line,  and  define 
the  state  space  of  Z  by  the  set  S  c  N,  the  set  of  natural 
numbers.  We  specifically  assume  that  Z  has  a  finite  state 


space  S  =  {1, . . . ,  K)  with  K  >2.  Let  R(t)  be  defined  as 
the  degradation  rate  of  the  system  at  time  t  e  R+  and  de¬ 
fine  a  positive  function  r  :  S  — >•  (0,  oo).  The  properties  of 
the  function  /■(•)  are  dictated  by  the  type  of  system  under 
consideration  and  its  surrounding  environment.  Since  the 
degradation  rate  of  the  system  is  explicitly  dependent  on  the 
environment  process,  the  degradation  rate  assumes  values 
in  the  space  V  =  {r(  1), . . . ,  r(K)}. 

Next,  we  describe  the  stochastic  evolution  of  the  sys¬ 
tem.  If  Z(t)  =  i  e  S,  then  (i)  R(t)  :=  r(Z(t))  =  r(i)  e  T>; 
and  (ii)  the  environment  transitions  from  state  i  e  S  to 
state  j  e  S  at  time  t  according  to  a  Markov  transition  func¬ 
tion  P(t)  :=  \pij(t)\  where  pUj(t)  :=  P{Z(t)  =j\Z(0)  =  /}. 
The  environment  process  Z  is  assumed  to  be  a  temporally 
homogeneous,  finite-state  Markov  process  so  that  pLJ(t) 
does  not  depend  on  t  for  all  /',  /  e  S.  Denote  by  X(t), 
the  cumulative  degradation  of  the  single-unit  system  up  to 
time  t  e  R+.  The  continuous- time,  continuous-state  degra¬ 
dation  process  [X(t)  :  t  >  0}  assumes  values  on  the  non¬ 
negative  real  line.  Moreover,  this  monotonically  increasing 
process  is  a  continuous,  additive  functional  of  Z,  and  thus, 
(X,  Z)  constitutes  a  special  case  of  a  Markov  additive  pro¬ 
cess  ((,’inlar,  1977).  The  main  contribution  of  Kharoufeh 
(2003)  was  a  closed-form  expression  for  the  cumulative 
distribution  function  (cdf)  and  moments  of  the  system 
lifetime  via  an  analysis  of  the  bivariate  Markov  process 
{(X(t),  Z(t))  :  /  >  ()}.  Next,  we  briefly  review  those  results 
without  proof.  Lor  specific  details  and  proofs,  the  reader  is 
referred  to  Kharoufeh  (2003). 

The  cumulative  degradation  of  the  single-unit  system 
up  to  time  /  e  1,  is  defined  by  the  cumulative  stochastic 
process: 

X(t)=  [  r(Z(u))du,  (1) 

Jo 

when  X(0)  =  0.  The  system  fails  as  soon  as  the  magnitude 
of  its  accumulated  degradation  exceeds  a  fixed  threshold 
value  x  (i.e.,  a  soft  failure).  The  lifetime  of  the  system  is 
given  by  the  random  variable: 

T(x)  =  in f{/  :  X(t)  >  x},  (2) 

or  the  first  time  the  X  process  exceeds  the  threshold  value 
x.  Define  the  probability  distributions: 

Vij(x,  t )  =  P{X(t)  <  x,  Z(0  =./|Z(0)  =  /},  ij  e  S , 

(3) 

where  F,j(x,  t)  is  the  joint  probability  that,  at  time  /,  the 
degradation  of  the  system  has  not  exceeded  the  value  x 
and  the  environment  process  is  in  state  j  e  S  given  that  the 
environment  was  initially  in  state  i  e  S.  The  distribution 
matrix  of  X{t)  is  defined  as: 

=  [Vij(x,t)\.  (4) 

Due  to  the  dual  relationship  of  Equation  (2),  the  cdf  of 
7’(x),  the  unit’s  full  lifetime,  is  obtained  as: 

F(x,  t)  \=  P{T(x)  </}  =  !  —  «V(x,  t)e,  (5) 
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where  a:  =  [a,]  with  a,:  =  P{Z( 0)  =  /}  and  e  is  a  K- 
dimensional  column  vector  of  ones. 

Of  particular  importance  to  the  problem  of  lifetime  prog¬ 
nosis  is  evaluation  of  the  unit’s  residual  lifetime  distribution 
given  that  it  has  not  failed  by  time  £o  >  0-  We  may  write  this 
distribution  directly  as: 

S(x,  t  |  ft)  :=  P{T{x )  >  t  +  £o  I  T(x)  >  &}, 

=  (1  -  F{x,  t  +  |o))/(l  -  F(x,  $o)), 

=  aV(x,  t  +  £0)e/(aV(x,  §o)e).  (6) 


By  integrating  the  tail,  the  nth  moment  of  the  unit  lifetime, 
for  n  >  1 ,  is  given  by: 


m(n)(x) 


E[Tn{x )]  =  / 


ntn~lP{T(x)  >  t}dt 
ntn~l(a\(x,  t)e)dt. 


(7) 


and  the  nth  moment  of  the  unit’s  residual  lifetime,  for  n  >  1 , 
is  given  by: 

m(n)(x  |  £0)  :=  E[T\x)  \  T(x)  >  $„], 

=  f  ntn~lS(x,  t  |  £o)dt, 

J  R+ 


aV(i 


i  r 

c,  £o)e  J%0 


nt"  1(aV(x,  t)e)dt.  (8) 


Equations  (6)  and  (8)  indicate  that  the  residual  lifetime  dis¬ 
tribution  and  moments  may  be  computed  via  the  full  life¬ 
time  distribution.  For  this  reason,  we  need  only  concern 
ourselves  with  the  numerical  evaluation  of  Equation  (5). 
Kharoufeh  (2003)  showed  that  the  matrix,  V(x,  t),  satisfies 
the  partial  differential  equation: 


Vt(x,  t)  +  Yx(x,  ORd  =  V(x,  t) Q,  (9) 


where  Yt(x,  t )  and  Vv(x,  / )  denote  the  partial  derivatives 
of  V(x,  i )  with  respect  to  t  and  x,  respectively,  = 
diag(r(l), . . . ,  r(K)),  and  Q  is  the  infinitesimal  generator 
matrix  for  {Z(t)  :  t  >  0}.  The  solution  to  the  PDE,  obtained 
in  the  two-dimensional  transform  space,  is  given  by: 

V*(«,  s)  =  (wRd  +  M  —  Q)_1,  Re(«)  >  0,  Refs)  >  0, 

(10) 

where  V*(«,  5)  is  obtained  by  the  two  matrix  transforms. 

V*(x,  s)  =  f  e~stY(x,  t)dt,  (11) 

Jr+ 

the  Laplace  transform  of  V(x,  t)  with  respect  to  t  and 

V*(u,  s)  =  f  e~uxY*(dx,  s),  (12) 

J  R+ 

the  Laplace-Stieltjes  Transform  (LST)  of  V*(x,  5)  with  re¬ 
spect  to  x.  Furthermore,  we  define  the  following  trans¬ 
forms: 

F*(x,s)=  f  e~s> F(x,  t)dt , 

J  IRq. 


F*(u,s)=  /  e~uxF*(dx,s). 

J  K+ 

The  full  lifetime  distribution  of  a  single-unit  system  in  a 
homogeneous  Markovian  environment  is: 

F*(u,  5)  =  —  aV*(w,  x)e,  Re(w)  >  0,  Refs)  >  0. 

(13) 

The  two-dimensional  transform  of  Equation  (13)  can  be  in¬ 
verted  numerically  using  inversion  algorithms  such  as  those 
due  to  Moorthy  (1995)  and  Abate  et  al.  (1998). 

Elsing  Equation  (13),  it  was  further  shown  by  Kharoufeh 
(2003)  that  the  (LST)  of  the  nth  moment  of  the  unit  lifetime, 
for  n  >  1 ,  is  given  by: 

hfn>(u)  =  n\a(uRu  —  Q)~"e.  (14) 

The  transform  of  Equation  (14)  can  be  inverted  numerically 
using  the  one-dimensional  inversion  algorithm  of  Abate 
and  Whitt  (1995).  Unfortunately,  a  closed-form  expression 
for  the  LST  of  the  nth  residual  lifetime  moment  (n  >  2)  does 
not  appear  to  exist.  Elowever,  the  mean  residual  lifetime  may 
be  obtained  by  evaluating  (numerically)  the  integral: 

w(I)(x  |  §0)  =  -TT- — —  /  («V(x,  t)e)dt.  (15) 
aV(x,  ^0)e  Jf(1 

The  analytical  models  of  this  section  provide  a  viable 
means  for  residual  lifetime  prognosis.  In  the  sections  that 
follow,  we  extend  the  basic  models  and  apply  them  to  two 
distinct  scenarios  in  which  real  sensor  data  may  be  used  to 
obtain  the  pertinent  measures.  Both  models  require  statis¬ 
tical  estimation  of  the  parameters  of  the  governing  Markov 
chain. 


3.  Model  I:  Observable  environment 

In  this  section,  it  is  assumed  that  sensor  data  provide  infor¬ 
mation  regarding  the  current  state  of  the  ambient  environ¬ 
ment  in  which  the  single-unit  system  resides  while  the  degra¬ 
dation  status  of  the  unit  is  unobservable.  Similar  models, 
developed  for  engineering  and  medical  applications,  can  be 
found  in  the  current  lifetime  analysis  literature  (Whitmore 
et  al.,  1998;  Lee  et  al.,  2000).  Our  environment  process  is 
analogous  to  a  “marker”  process  whereas  the  degradation 
process  corresponds  to  a  “latent”  process  that  is  assumed 
to  be  unobservable.  In  those  works,  the  marker  and  latent 
processes  were  assumed  to  form  a  bivariate  Brownian  mo¬ 
tion  process.  In  this  paper,  we  impose  the  assumption  that 
the  marker  process  is  a  Markov  chain  that  directly  influ¬ 
ences  the  latent  (degradation)  process;  however,  we  make 
no  assumptions  regarding  the  probability  law  of  the  latter. 

Under  its  prevailing  operating  conditions,  the  system 
wears  continuously  (and  additively)  until  its  cumulative 
degradation  exceeds  a  fixed,  deterministic  value  x  at  which 
time  the  unit  is  said  to  have  failed.  Optimally-located  sen¬ 
sors  provide  real-time  data  regarding  the  conditions  of  the 
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environment  (e.g.,  temperature,  pressure,  humidity,  etc.); 
however,  the  cumulative  amount  of  degradation  cannot  be 
discerned  from  these  readings.  It  is  further  assumed  that 
the  environment  state  space  may  be  partitioned  into  K  dis¬ 
tinct  states  as  S  =  {1, 2, ... ,  K},  K  <  oo.  Associated  with 
each  distinct  environment  state  is  a  known  degradation  rate, 
namely  r  (/'),  /'  e  S  so  that  when  the  environment  is  in  state 
j  e  .S’,  the  system  degrades  at  rate  r(j).  Knowledge  of  the 
real-valued  function  r  is  assumed  to  be  available  through 
physical  properties  of  the  unit  under  consideration.  Using 
this  information,  it  is  possible  to  directly  apply  the  results  of 
Section  2  to  characterize  the  full  and  residual  lifetime  distri¬ 
butions.  However,  in  order  to  do  so,  we  require  a  surrogate 
(estimated)  stochastic  process  for  the  true,  underlying  envi¬ 
ronment  process. 

Equation  (13)  requires  specification  of  three  key  com¬ 
ponents:  (i)  the  initial  distribution  vector  (a)  for  the  en¬ 
vironment  process:  (ii)  the  diagonal  matrix  of  degradation 
rates  (Rd);  and  (iii)  the  infinitesimal  generator  matrix  (Q) 
for  the  environment  process.  Since  the  environment  process 
is  observable  and  the  degradation  rates  are  assumed  to  be 
known  for  each  environmental  state,  it  is  necessary  only  to 
approximate  the  matrix  Q.  We  do  so  by  employing  stan¬ 
dard  statistical  inference  techniques  for  Markov  processes 
(Basawa  and  Rao,  1980). 

We  assume  that  the  environment  process  \Z(i )  :/>()}  is 
observable  (and  homogeneous)  over  the  time  interval  [0,  J], 
The  system  has  not  failed  by  time  T,  and  began  its  lifetime 
in  perfect  working  order  at  time  0  (i.e.,  X(0)  =  0).  The  en¬ 
vironment  is  continuously  observed  up  to  time  T  so  that,  at 
each  transition  epoch,  we  record  the  current  and  subsequent 
states  of  the  random  environment.  Let  q(i,j)  denote  the  rate 
at  which  the  environment  transitions  from  state  i  e  S  to 
state  j  e  S,j  /  i.  If  the  environment  process  {Z(t)  :  t  >  0} 
is  a  continuous-time  Markov  chain,  then  the  holding  time 
in  state  i,  given  that  the  subsequent  state  is  j  /  i,  is  expo¬ 
nentially  distributed  with  rate  parameter  q(i.j).  Let  NT(L  j) 
denote  the  random,  integer  number  of  transitions  of  the 
process  from  i  to  j  in  a  time  interval  of  length  T.  Moreover, 
let  II j(i)  denote  the  total  holding  time  in  state  i  e  S  during 
[0,  T ].  It  can  be  shown  (Basawa  and  Rao,  1980)  that: 

qCj)  =  lim  ■  (16) 

T->  oo  HT(l) 

Therefore,  for  sufficiently  large  T,  we  may  approximate 
/  i  by: 

q  (/,  j)  qT(i,  j )  =  •  ( 1 7) 

Ht{i) 

The  diagonal  elements  of  the  generator  matrix  are  obtained 
as: 

qT{i,  /)  =  -  faiUj),  i  £  S.  (18) 

j¥=i 

We  pause  here  to  note  that,  for  a  fixed  observation  in¬ 
terval,  it  is  possible  to  alternatively  estimate  the  generator 


matrix  by  observing  k  independent  sample  paths  of  the  Z 
process.  Define  N(k)(i,  j)  as  the  total  number  of  transitions 
from  state  i  to  j  over  all  k  trials  and  define  H(k)(i)  as  the  to¬ 
tal  holding  time  in  state  i  over  all  k  trials.  The  off-diagonal 
elements  of  the  generator  matrix  are  given  by: 


q('J)  = 


lim 

k—>oo 


H(kHi)  ’ 


(19) 


so  that  for  k  sufficiently  large,  the  approximation  of  q(i,j), 
j  /  i  is: 


q(i,j)  w  q{  \i,j )  = 

We  obtain  the  diagonal  elements  by: 

$(*>(/,  0  =  -£y*} (i,j)  ieS. 

Mi 


(20) 

(21) 


It  has  been  shown  in  Basawa  and  Rao  (1980)  that: 

,.  .,  ,.  NT{Uj)  v  N«\Uj) 

q(i,j)  =  Inn  „  —  =  Inn  .  . 

T-*  oo  HT(l)  oo  7/A)p) 


(22) 


In  model  I,  we  observe  a  single  sample  path  over  a  suffi¬ 
ciently  long  period  using  Equations  (17)  and  (18),  whereas 
we  utilize  Equations  (20)  and  (21)  in  model  II.  Let  Q  j  := 
[c/tO-J)]  and  let  V(x,  / )  denote  the  matrix  solution  to  Equa¬ 
tion  (9)  when  the  generator  matrix  QT  is  used  as  a  surro¬ 
gate  for  Q.  The  estimated  matrix  distribution,  V(x,  t),  is 
obtained  through  the  double  Laplace  transform  inversion: 


V(x,  t)  =  jC.-1{u~\uRb  +  s1  -  Q rr‘}, 

where  Re(t<)  >  0,  Refs)  >  0,  and  C  1  denotes  the  (double) 
inverse  Laplace  operator. 

The  resulting  estimates  for  the  full  and  residual  lifetime 
distributions  are,  respectively. 

F(x,  t)=  1  -  txV(x,  ?)e,  (23) 


and 


S(x,  t  |  ^0)  = 


<xV(x,  t  +  g0)e 
ccV(x,  ^0)e 


(24) 


We  illustrate  the  estimation  procedure  via  numerical  exam¬ 
ples  in  Section  5.  In  model  II,  we  assume  that  only  degra¬ 
dation  measures  are  observable  for  the  single-unit  system. 


4.  Model  II:  Observable  degradation 

The  ultimate  purpose  of  a  prognostic  system  is  to  assess 
the  current  “health”  of  the  unit  and  to  make  inferences  re¬ 
garding  the  remaining  useful  lifetime  of  the  unit.  Model  II 
provides  a  viable  first  step  towards  an  analytical  approach 
to  address  this  problem  when  the  component’s  degradation 
is  governed  by  a  stationary,  Markov  environment.  Assume, 
at  time  t  >  0,  that  the  cumulative  degradation  of  the  single¬ 
unit  system,  X{t),  is  observable.  In  contrast  to  model  I,  the 
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environment,  its  number  of  distinct  states,  and  its  degrada¬ 
tion  rates  are  unknown. 

The  basis  of  our  procedure  can  be  described  as  follows. 
We  approximate  the  true  degradation  path  using  simple, 
piecewise  linear  functions  that  pass  through  a  finite  num¬ 
ber  of  observations.  The  slope  of  each  line  segment  approx¬ 
imates  the  degradation  rate  over  a  given  time  interval.  Once 
the  slopes  have  been  collected,  we  perform  -means  cluster 
analysis  wherein  the  centroid  of  each  cluster  corresponds 
to  the  mean  degradation  rate  for  a  distinct  state  of  the  envi¬ 
ronment.  However,  selection  of  the  integer  K  is  non-trivial. 
We  describe  a  means  by  which  to  estimate  this  value  in  the 
next  subsection. 


4.1.  Estimating  the  number  of  states 

In  the  absence  of  environmental  state  observations,  it  is 
difficult  to  discern  the  total  number  of  distinct  states  that 
might  be  encountered  by  the  unit.  In  some  applications,  the 
number  of  states  might  be  obvious.  For  example,  in  a  man¬ 
ufacturing  context,  if  a  cutting  tool  operates  at  only  three 
distinct  cutting  speeds  (e.g.,  slow,  moderate,  fast),  and  the 
wear  rate  of  the  cutting  tool  is  known  for  each  speed,  then 
three  distinct  “environment”  states  exist.  In  other  applica¬ 
tions,  however,  it  may  be  necessary  to  rely  on  the  knowledge 
and  experience  of  subject  matter  experts  to  determine  an 
appropriate  number  of  environment  states.  In  this  subsec¬ 
tion,  we  address  scenarios  in  which  the  appropriate  integer 
number  of  states  is  not  obvious,  and  the  engineer  has  at 
his/her  disposal  only  the  observed  degradation  of  the  unit. 
The  problem  is  that  of  selecting  an  appropriate  partition  of 
the  set  of  all  observed  degradation  rates  over  some  finite 
observation  period. 

Fraley  and  Raftery  (1998)  note  that  the  method  of  K- 
means  clustering  is  the  most  widely  used  non-hierarchical 
method  for  partitioning  a  set  of  real  observations  into  K 
distinct  clusters  if  the  value  K  is  known  a  priori.  We  now 
describe  a  method  to  estimate  the  appropriate  value  of  K 
when  a  total  of  N  degradation  rates  have  been  observed. 
The  technique  performs  a  comparison  of  K  means  via  a 
one-way  analysis  of  variance  by  computing  an  appropriate 
/•'-ratio  that  depends  explicitly  on  K. 

Denote  by  y,  the  mean  of  cluster  i,  i  =  1,  2, . . . ,  K,  and 
let  77,  denote  the  number  of  observations  in  the  7th  cluster- 
such  that: 

K 

*  =  £»<• 

7=1 

The  overall  sample  average  of  degradation  rates  is: 


Define  by  SSB,  the  sum  of  squares  between  clusters  which 
is: 

K 

SSB=Yjn,iyi-y)2.  (25) 

7=1 

Similarly,  the  sum  of  squares  within  clusters  (SS  W)  is: 

K  rii 

ssw  =  YJY.(yd-yrf'  w 

7=1  7=1 

where  ytJ  denotes  the  / th  degradation  rate  observed  in  the 
/th  cluster.  The  /-'-ratio  for  K  clusters,  2  <  K  <  N,  is: 

r  SSB/(K-  1) 

Fk  =  SSWUN-K )'  <27) 

Our  objective  is  to  find  the  “best”  value  of  K  via  Equa¬ 
tion  (27).  Calinski  and  Harabasz  (1974)  suggest  choosing 
K  such  that,  over  some  set  1C,  Fk  is  the  absolute  maxi¬ 
mum,  the  first  local  maximum,  or  a  point  at  which  the 
function  exhibits  a  comparatively  rapid  increase.  It  is  worth 
noting  that,  if  FK  is  monotonically  increasing  in  K,  then 
the  number  of  clusters  (states)  should  be  equal  to  the  total 
number  of  observations  (i.e.,  each  observation  constitutes 
a  cluster).  Let  K.  =  {2,  3, . . . ,  /}  where  J  is  a  positive  inte¬ 
ger.  For  the  sake  of  computational  expedience,  it  is  ideal 
to  choose  the  smallest  value  of  K  that  leads  to  an  appro¬ 
priate  representation  of  the  underlying  process.  As  a  first 
course  of  action,  we  apply  the  second  criterion  suggested 
in  Calinski  and  Harabasz  (1974),  and  choose  a  value  of  K 
that  corresponds  to  the  first  local  maximum  over  K.  given 
by: 

k  =  mm{K  elC:  Fk>  Fk+1).  (28) 

However,  in  the  case  where  FK  is  strictly  increasing  in  K,  we 
resort  to  the  remaining  two  criteria  to  estimate  the  small¬ 
est  possible  value.  We  illustrate  and  apply  the  estimation 
procedure  in  Section  5. 

4.2.  Description  of  the  estimation  procedure 

A  formal  description  of  the  full  estimation  procedure 
(model  II)  is  now  provided. 

Step  0.  Initialization. 

At  time  ?o  =  0,  observe  A(70)  =  x0. 

Step  1.  Observe  degradation  measures. 

Observe  the  degradation  at  times  ?i  <  •  •  •  <  7m, 
M  e  N  and  form  the  set  of  observation  times: 

7~  :={(/:/  =  0,  1, 2, ... ,  M).  (29) 

with  observations  X(tj),j  =  0,  1,  2, _ M.  It  is  as¬ 

sumed  that,  at  time  7m,  the  single-unit  system  has 
not  failed. 

Step  2.  Approximate  the  failure  path. 

After  observing  the  degradation  path  up  to  time 
epoch  7m,  approximate  the  true  failure  path  by 
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a  simple  piecewise-linear  approximation  that  con¬ 
nects  the  observed  degradation  measures  for  each 
element  of  T. 

Step  3.  Approximate  degradation  rates  via  finite  difference 
methods. 

For  each  observation  time  tj  G  T,  approximate  the 
potential  degradation  rates  of  the  process  by  the 
difference  equation: 


Yj  ■= 


Xjtj)  ~  X(tj- 1) 
h  —  0-i 


j  =  1,2, ,  M.  (30) 


5.1.  Model  I  examples 

For  each  j  e  S.  there  exists  a  known  degradation  rate  for  the 
single-unit  system,  namely  r(j).  The  environment  evolves 
over  time  and  is  observed  over  the  interval  [0,7],  We  pro¬ 
vide  three  test  cases,  namely  when  the  environment  process 
(Z)  assumes  values  on  a  state  space  with  |Sj  =  2,  5  and  10, 
respectively.  For  each  test  case,  the  elements  of  the  generator 
matrix  Q  were  drawn  from  a  continuous  uniform  popula¬ 
tion  on  (20,  40)  such  that,  for  j  /  i,  q(i,j)  ~  17(20,  40).  In 
the  usual  way,  the  diagonal  elements  of  the  generator  matrix 
are: 


We  assume  that  the  inter-observation  times  (fi  — 
tj- 1,  j  =  1,2,...,  M)  are  equal.  For  the  discrete 
sampling  interval,  [to.  f  v/],  gather  the  observed  rates 
in  a  set: 

F  =  {»■  :y  =  1, 2, . . . ,  M). 

Step  4.  Compute  the  K  distinct  degradation  rates. 

Select  an  appropriate  value  K  using  the  proce¬ 
dure  outlined  in  Section  4. 1  applied  to  observations 
in  the  set  F.  This  may  be  accomplished  by  using 
any  standard  statistical  software  package.  Define 
C  =  {Ci,  Ci, _ Cf.)  as  the  set  of  K  distinct  clus¬ 

ters  such  that  Cj  n  C)  =  0,y  /  and  /q  denotes  the 
centroid  of  cluster  C(,  i  =  1,2 , ...  ,K.  Each  yj  g  T 
is  therefore  assigned  to  exactly  one  cluster  in  C  such 
that  the  estimated  degradation  rate  of  environment 
state  j  g  S  is  f  (/)  =  /x7. 

Step  5.  Construct  new  degradation  path. 

Construct  a  new  piecewise  linear  degradation  path 
by  replacing  each  y7  g  V  by  f(j)  =  / ij .  (Note:  if 
r(j)  =  r(j  +  1)  for  some  j,  then  the  two  adjacent 
line  segments  are  replaced  by  a  single  line  segment 
with  slope  equal  to  r  (/).) 

Step  6.  Approximate  generator  matrix. 

Using  the  piecewise  linear  estimate  of  the  degra¬ 
dation  path  of  Step  5,  estimate  the  ^-dimensional 
generator  matrix  Q  by  Q,w  using  Equations  ( 1 7) 
and  (18). 

Remark  1.  Three  important  assumptions  are  employed 

here:  (i)  degradation  of  the  unit  is  perfectly  observable; 

(ii)  at  time  tM,  the  single-unit  system  has  not  failed;  and 

(iii)  the  degradation  paths  are  monotonically  increasing. 


5.  Numerical  examples 

This  section  provides  illustrative  examples  of  the  failure 
time  models  (I  and  II)  of  the  two  preceding  sections.  We 
validate  our  technique  by  comparing  the  full  and  residual 
lifetime  probability  values  with  simulated  values.  Moreover, 
we  test  model  II  by  comparing  results  with  empirical  ob¬ 
servations  of  unit  lifetimes. 


q(f  i)  =  -  ^q(f  j),  i  G  .S'. 

j¥=i 

Moreover,  the  true  generator  matrix  Q  is  estimated  by 
the  matrix  Q  j  with  T  =  100,  500,  5000  using  Equations 
(17)  and  (18).  The  associated  degradation  rates  (r(j),j  = 
1,2 , ...  ,K)  were  drawn  from  continuous  uniform  popula¬ 
tions  on  the  interval  (20,  80).  In  all  test  cases,  the  failure 
threshold  value  was  fixed  at  .v  =  1.0  units.  Full  and  resid¬ 
ual  probability  values  were  computed  via  Equations  (5)  and 
(6),  respectively,  using  a  variant  of  the  inversion  algorithm 
in  Moorthy  (1995). 

We  compare,  at  m  fixed  points,  r1?  r2, . . . ,  rm,  distribu¬ 
tion  function  values  generated  by  the  true  and  estimated 
processes.  Goodness-of-ht  tests  were  conducted  (at  the  0.05 
level)  using  the  Cramer- von  Mises  test  statistic  to  compare 
full  and  residual  lifetime  distributions.  The  appropriate  test 
statistic,  denoted  by  K2,  is  given  by: 

i  m 

K2  =  z£(Gi(tO-G2(tO)2,  (31) 

1  i=  l 

where  G\  denotes  the  distribution  function  obtained  from 
the  estimated  environment  process  and  (f  is  the  cdf  gen¬ 
erated  using  the  true  environment  process.  We  denote  the 
critical  value  by  k*.  If  ki  <  k*.  then  we  fail  to  reject  the  null 
hypothesis  that  G\  =  63.  Table  1  provides  a  summary  of 
the  Cramer-von  Mises  test  statistics  when  comparing  cu¬ 
mulative  probability  values  obtained  by  simulating  the  true 
process  with  generator  Q  with  those  of  numerical  Laplace 
transform  inversion  obtained  using  Qj.  Three  observation 
periods  were  considered  in  this  experiment.  For  the  residual 
life  tests,  we  fixed  =  E[T{x)]. 

It  is  noted  that,  in  all  18  experiments,  we  fail  to  reject 
the  null  hypothesis  that  the  estimated  and  true  distribution 
functions  are  equivalent.  Some  additional  remarks  are  war¬ 
ranted  regarding  these  results. 

1.  If  the  environment  process  can  be  partitioned  into  K 
distinct  and  observable  states,  our  approach  provides  a 
viable  approximation  procedure  that  does  not  require 
failure  time  observations;  it  requires  only  a  count  of  en¬ 
vironment  transitions. 

2.  In  the  case  where  the  environmental  state  holding 
times  are  non-exponentially  distributed,  the  procedure 
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Table  1.  Cramer-von  Mises  test  statistics  for  model  I  experiments  (k*  =  0.461,  a  =  0.05) 


Distribution 

Run  length 

<N 

II 

■<-> 

II 

K  =  10 

F(x,  t) 

T  =  100 

7.013  x  10“3 

4.233  x  10-2 

9  All  x  10“3 

OO 

'xT 

II 

II 

O 

o 

1.649  x  10“3 

7.895  x  10-3 

2.130  x  10-3 

T  =  5000 

9.941  x  10“5 

9.019  x  10“5 

1.986  x  10“6 

1  -  S(x,  r|$0) 

T=  100 

5.511  x  10-3 

3.821  x  10-2 

7.455  x  10-3 

(m  =  50) 

II 

o 

o 

1.223  x  10“3 

6.553  x  10“3 

2.113  x  10-3 

T  =  5000 

5.544  x  10-5 

8.092  x  10-5 

2.576  x  IQ”6 

remains  feasible  since  the  holding  time  distributions  may 
be  approximated  by  phase-type  distributions  (e.g.,  gen¬ 
eralized  Erlang  distributions)  that  retain  the  memoryless 
property.  We  shall  elaborate  further  on  this  point  in  Sec¬ 
tion  6. 


5.2.  Model  II  examples 

In  this  subsection,  we  first  illustrate  model  II  using  simu¬ 
lated,  linear  degradation  paths  before  testing  the  procedure 
on  a  set  of  non-linear  degradation  data.  For  the  simulation 
experiment,  sample  paths  were  generated  via  a  known  en¬ 
vironment  process  with  K  distinct  states  and  known  degra¬ 
dation  rates  {r(l),  r(2),  . . . ,  r(K)}.  Initially,  we  fix  K  =  3. 
Five  hundred  degradation  sample  paths,  similar  to  the  five 
shown  in  Fig.  1,  were  generated  for  the  experiment. 

Setting  to  =  0,  we  observe  the  level  of  degradation  at  M 
equally  spaced  points  in  time  such  that  tM  =  20.  That  is, 
the  degradation  process  is  observed  up  to  time  20.  Figure  2 
depicts  a  sample  of  observations  ( M  =  10)  obtained  from 
the  five  degradation  sample  paths  of  Fig.  1  that  were  used 
to  approximate  a  piecewise,  linear  degradation  sample  path 
by  the  procedure  of  Section  4. 

Next,  we  utilize  the  approximation  to  obtain  the  full  and 
residual  lifetime  distributions.  In  each  case,  we  observe  the 
system  up  to  time  Im  =  20  except  that  we  vary  the  inter¬ 


observation  time.  In  particular,  we  consider  M  =  20,  100 
and  200  observations  on  [0,  20],  respectively.  The  assumed 
unit  failure  threshold  was  fixed  at  x  =  1 52.6  units.  The  same 
500  simulated  degradation  sample  paths  were  used  for  each 
value  of  M  to  estimate  a  new  generator  matrix  and  new 
degradation  rates  at  times  4.0,  8.0  and  12.0  respectively. 

As  before,  we  compare  the  estimated  versus  “true”  full 
lifetime  cdfs  via  the  Cramer-von  Mises  goodness-of-ht  test. 
Additionally,  we  compare  the  estimated  versus  true  resid¬ 
ual  lifetime  distribution.  It  should  be  noted  that,  in  order  to 
test  the  null  hypothesis,  we  use  the  cdf  of  the  residual  life¬ 
time  distribution  (1  -  S(x,  t|£0))-  The  summarized  results 
of  Table  2  indicate  that  none  of  the  tests  are  significant  at 
the  0.05  level;  thus,  we  are  able  to  adequately  estimate  the 
full  and  residual  lifetime  distributions  for  these  simulated 
degradation  paths. 

To  illustrate  our  method  to  select  an  appropriate  number 
of  states  (A),  we  simulated  multiple  degradation  sample 
paths  of  environment  processes  having  five  and  10  states, 
respectively.  The  /-'-ratios  for  each  case  are  displayed  in 
Table  3.  Applying  the  first  local  maximum  rule,  the  /-'-ratios 
suggest  K  =  4  for  the  five-state  environment  and  K  =  8  for 
the  10-state  environment. 

Consequently,  we  compared  the  full  and  residual  lifetime 
distributions  using  K  and  K  +  1  states  in  each  scenario. 
Table  4  indicates  that  we  fail  to  reject  the  null  hypotheses 


Fig.  1.  A  sample  of  five  linear  degradation  paths. 


Fig.  2.  Piecewise-linear  approximation  of  the  degradation  paths. 


540 


Kharoufeh  and  Cox 


Table  2.  Cramer- von  Mises  test  statistics  for  simulated  data  (k*  = 
0.461  at  a  =  0.05) 


M 

Interval 

F(x,  t ) 

Ho 

1  -  S(x,  t  |  So) 

20 

[0. 0,4.0] 

0.0157 

4.0 

0.0189 

[0.0, 8.0] 

0.0210 

8.0 

0.1213 

[0.0,12.0] 

0.0247 

12.0 

0.1633 

too 

[0. 0,4.0] 

0.0095 

4.0 

0.0197 

[0.0, 8.0] 

0.0088 

8.0 

0.0476 

[0.0,12.0] 

0.0087 

12.0 

0.0639 

200 

[0. 0,4.0] 

0.0053 

4.0 

0.0149 

[0.0, 8.0] 

0.0051 

8.0 

0.0336 

[0.0,12.0] 

0.0052 

12.0 

0.0266 

Table  4.  Cramer-von  Mises  test  statistics  with  K  states  (k*  = 
0.461,  a  =  0.05,  m  =  100) 


K 

K 

F(  x ,  t ) 

Ho 

l-S(xJ\So) 

5 

4 

2.3056  x  10“2 

17.5311 

1.0179  x  lO”1 

5 

5 

3.2317  x  10-2 

17.5311 

1.0733  x  lO”1 

10 

8 

3.5704  x  10-2 

13.4249 

1.6783  x  lO”1 

10 

9 

3.5397  x  10-2 

13.4249 

1.4709  x  lO”1 

whenever  the  crack  length  first  exceeds  a  critical  value  of 
x  =  45  mm.  We  estimate  the  off-diagonal  generator  matrix 
values  by: 


that  the  full  and  residual  lifetime  distributions  are  equiva¬ 
lent  at  the  0.05  level. 

Finally,  we  illustrate  model  II  using  real  degradation  data 
which  correspond  to  the  propagation  of  a  fatigue  crack  in 
metallic  materials.  The  data  were  originally  analyzed  by 
Virkler  et  a!.  (1979)  who  measured  crack  length  propaga¬ 
tion  in  68  test  specimens  of  a  2024-T3  aluminum  alloy.  The 
crack  length  of  each  sample  was  measured  over  time  (in 
number  of  load  cycles).  Figure  3  depicts  a  representative 
sample  of  five  degradation  paths  from  the  original  data  set. 
For  a  more  thorough  treatment  of  fatigue  crack  dynam¬ 
ics,  the  reader  is  referred  to  Ray  and  Tangirala  (1997)  and 
references  contained  therein. 

To  illustrate  our  procedure,  let  X(t)  denote  the  length 
of  the  crack  at  time  t  and  assume  the  rate  at  which  the 
crack  grows  is  subject  to  its  random  environment  (ap¬ 
plied  stress,  ambient  conditions,  and  other  factors).  It  is 
assumed  that  these  environmental  factors  can  be  charac¬ 
terized  by  {Z(t)  :  t  >  0},  a  temporally  homogeneous  (time¬ 
stationary)  Markov  process  on  a  finite-state  sample  space. 
Although  the  underlying  data  may,  in  fact,  exhibit  non¬ 
stationary  behavior,  we  shall  observe  all  68  sample  paths 
to  estimate  the  generator  matrix  (Q)  for  the  Z  process.  In 
the  experiment  that  follows,  the  component  is  said  to  fail 


Table  3.  F-ratio  values  for  five-  and  10-state  environments  ( x  10s) 


K 

Five-state 

10-state 

2 

1.654973 

0.285050 

3 

3.090  848 

0.400  770 

4 

4.067  369 

0.421  102 

5 

3.820  857 

0.660245 

6 

3.249  939 

0.698  204 

7 

6.678  639 

0.955  170 

8 

7.278  563 

1.144632 

9 

10.76289 

1.047  932 

10 

7.383  198 

0.697  347 

11 

5.876  550 

1.306  396 

12 

9.920  737 

1.805007 

13 

8.558  721 

0.857228 

14 

15.161074 

0.858  865 

15 

12.455  254 

1.386817 

q(Uj) 


.(68)  =  N^KiJ) 
Ht6S)(i)  ’ 


(32) 


as  defined  in  Equation  (20).  It  should  be  noted,  however, 
that  the  empirical  cumulative  probabilities  for  the  data  were 
computed  with  a  sample  of  only  68  observations.  Never¬ 
theless,  we  take  this  empirical  distribution  to  be  the  “true” 
distribution. 

The  /-'-ratio  test  for  selecting  K  was  applied  to  the  data 
by  observing  55,  77  and  91%  of  the  lifetime,  respectively. 
We  select  two  clusters  at  55%  of  the  lifetime,  nine  clusters 
at  77%  of  the  lifetime  (due  to  the  sharp  increase),  and  12 
clusters  at  91%  of  the  lifetime  as  seen  in  Table  5. 

Estimates  of  the  generator  matrix  and  degradation  rates 
were  constructed  for  K  =  2,  3,  4.  9  and  12  when  observing 
up  to  91%  of  the  lifetime  of  each  specimen.  The  resulting 
cdfs  (full  and  residual  lifetime)  for  each  case  were  compared 
to  the  empirical  distribution  at  a  fixed  number  of  points 
(m  =  65)  using  the  Cramer-von  Mises  test.  The  results  of 
the  numerical  tests  are  summarized  in  Table  6. 

The  results  for  the  full  lifetime  distributions  suggest  that 
for  K  =  2,  3  and  4,  we  fail  to  reject  the  null  hypothesis  that 


Fig.  3.  The  propagation  of  the  fatigue  crack  length  for  five  test 
specimens. 


Stochastic  models  for  degradation-based  reliability 
Table  5.  F- ratio  comparison  (x  104) 


Percentage  of  lifetime 


K 

55% 

77% 

91%, 

2 

0.845  912 

1.716694 

1.934667 

3 

0.770  946 

2.045  060 

2.448  228 

4 

0.700  849 

2.267  842 

2.866000 

5 

0.637918 

2.412888 

3.248  069 

6 

1.789474 

2.506080 

3.620  849 

7 

1.984733 

2.627177 

3.934  574 

8 

2.230  180 

2.658  314 

4.208  318 

9 

2.570793 

2.669  967 

4.374201 

10 

2.785035 

4.566412 

4.448  147 

11 

3.017  830 

4.974  103 

4.744162 

12 

2.993  362 

5.272  737 

4.830633 

13 

3.714449 

5.683212 

4.796  518 

14 

4.167025 

5.850034 

4.800008 

15 

4.292031 

6.608  150 

4.745  232 

the  distributions  are  equal.  For  the  residual  lifetime  dis¬ 
tributions  with  K  =  3  and  4,  we  fail  to  reject  the  null  hy¬ 
pothesis  that  the  distributions  are  equal.  The  /-'-ratio  test  of 
Table  5  indicates  an  estimate  of  A'  =  12.  However,  using  an 
estimated  12-state  generator  matrix,  the  resulting  full  and 
residual  lifetime  cdfs  fail  to  pass  the  goodness-of-ht  tests. 

Although  it  is  apparent  that  there  may  exist  greater  than 
three  or  four  environmental  states,  we  surmise  that  the  12- 
state  estimate  is  poor  for  a  few  reasons.  First,  to  adequately 
estimate  a  12-state  generator  matrix,  a  fairly  large  number 
of  environmental  transitions  must  occur  during  a  Unite  ob¬ 
servation  period.  It  is  obvious,  by  inspection  of  Fig.  3,  that 
few  significant  transitions  take  place  early  on.  Second,  it 
is  likely  that  transitions  occur  primarily  between  a  subset 
of  the  overall  environment  sample  space  at  different  times. 
That  is,  some  state  transitions  appear  to  dominate  early 
in  the  specimen  lifetimes  whereas  others  dominate  as  the 
specimens  age  (i.e.,  as  the  crack  grows).  Consequently,  we 
surmise  that  an  age-  or  degradation-state-dependent  model 
may  be  more  appropriate  for  this  data  set,  as  it  appears 
to  violate  the  time-stationarity  assumptions  of  our  model. 
Flowever,  it  is  important  to  note  that  the  estimation  proce¬ 
dure  indicated  that  only  a  few  states  could  be  used  within 
our  framework  to  provide  adequate  cumulative  probability 
estimates.  We  surmise  that  there  may  exist  a  minimal  rep- 

Table  6.  Cramer-von  Mises  test  statistic  (k*  =  0.461  at  a  =  0.05, 
m  =  65) 


A  p(x,t)  &  l-§(x,t\$o) 


2 

1.398  488  x  10-1 

2.534 

6.124815  x  lO”1 

3 

8.671  149  x  10“2 

2.534 

1.395  285  x  10_1 

4 

4.525  246  x  lO”1 

2.534 

4.383  202  x  10”1 

9 

7.657227  x  10”1 

2.534 

2.235  864 

12 

8.723  972  x  10”1 

2.534 

1.667  689 

541 


Table  7.  Mean  full  and  residual  lifetimes  ( x  105  load  cycles),  £o  = 
actual  m(X){x) 


x  ( mm ) 

if1) 

(x) 

if1) 

(x  \  Ho) 

Actual 

Model 

Actual 

Model 

10 

0.318  574 

0.325  974 

0.349  943 

0.347018 

15 

1.183  405 

1.150109 

1.268  688 

1.293  640 

20 

1.638  693 

1.530  864 

1.738  733 

1.771  963 

25 

1.938014 

1.800  568 

2.042441 

2.113  597 

30 

2.158  575 

2.021  006 

2.286  812 

2.368  369 

35 

2.326253 

2.116553 

2.458462 

2.508  772 

41 

2.464467 

2.298  747 

2.611653 

2.633  815 

resentation  for  the  underlying  process  that  includes  only 
the  “most  important”  states.  In  the  current  framework,  our 
model  does  not  dynamically  dictate  which  states  are  most 
important. 

Although  the  distribution  results  were  comparatively  dis¬ 
appointing  for  the  empirical  data,  we  conducted  an  experi¬ 
ment  to  compare  the  actual  mean  lifetime  and  mean  resid¬ 
ual  lifetime  with  the  results  of  our  models  for  various  failure 
threshold  levels.  We  computed  the  mean  residual  lifetimes 
by  conditioning  on  the  event  that  a  specimen  survives  be¬ 
yond  the  observed,  unconditional  mean  lifetime.  For  each 
threshold  value,  a  new  generator  matrix  and  degradation 
rates  were  estimated  using  K  =  3  environmental  states.  The 
mean  lifetimes  over  the  68  degradation  paths  were  com¬ 
pared  to  those  computed  by  Equation  (14).  The  mean  resid¬ 
ual  lifetimes  were  computed  by  Equation  (15).  The  initial 
measured  crack  length  was  approximately  9.0  mm  (i.e., 
A(0)  9.00  mm).  Table  7  summarizes  the  results  of  this 

experiment. 

The  results  of  this  experiment  indicate  that  the  procedure 
adequately  estimates  the  mean  lifetimes  and  the  mean  resid¬ 
ual  lifetimes.  Although  not  as  informative  as  the  residual 
life  distribution,  the  lower  moment  approximations  may  be 
useful  for  constructing  surrogate  parametric  distributions. 

6.  Conclusions 

In  this  paper,  we  have  presented  a  novel  approach  for  the 
estimation  of  full  and  residual  lifetime  distributions  (and 
moments)  for  single-unit  systems  subject  to  a  Markov  en¬ 
vironment.  The  approach  provides  a  necessary  first  step 
toward  a  formal,  analytical  technique  for  remaining  life¬ 
time  prognosis  via  real  degradation  measures.  Knowledge 
of  the  residual  lifetime  distribution  is  especially  useful  for 
prescribing  policies  that  may  reduce  the  risk  of  catastrophic 
failures  while  eliminating  superfluous  preventive  mainte¬ 
nance  activities. 

Our  approach  is  novel  in  that  we  combine  analytical 
stochastic  modelling  techniques  with  real  sensor  data  in  two 
distinct  cases.  First,  we  considered  the  case  when  the  sensor 
data  provide  real-time  information  regarding  the  current 
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state  of  the  environment  in  which  the  single-unit  system  re¬ 
sides.  Second,  we  considered  the  case  in  which  the  sensors 
provide  the  current  level  of  degradation  of  the  component. 
In  the  first  case,  our  numerical  results  show  great  promise 
for  a  fully-automated  technique  for  the  estimation  of  life¬ 
time  distributions  under  our  problem  assumptions.  In  the 
second  case,  we  demonstrated  that  our  technique  could  be 
used  to  compute  the  distributions  given  nothing  more  than 
degradation  measures  (assumed  to  be  perfectly  observable). 
Although  our  procedure  for  estimating  residual  lifetime  dis¬ 
tributions  performed  only  moderately  well  on  the  empirical 
data  of  Virkler  et  al.  (1979),  the  comparison  of  mean  resid¬ 
ual  lifetimes  is  very  promising. 

The  techniques  presented  in  this  research  have  raised  a 
few  important  questions  that  require  further  study.  First,  it 
may  be  possible  that  the  time  between  environment  tran¬ 
sitions  is  distinctly  non-exponential.  In  such  a  case,  the 
engineer  should  resort  to  phase-type  approximations  for 
the  state  holding  time  distributions  as  outlined  by  Altiok 
(1985).  This  approach  has  the  benefit  of  retaining  the 
Markov  property  of  the  environmental  process  so  that  the 
main  results  of  this  paper  (i.e.,  Equations  (5)  and  (6))  may 
still  be  applied.  The  disadvantage  of  this  approach  is  that 
it  requires  an  expansion  of  the  environment  state  space. 
Second,  there  is  the  issue  of  non-stationarity  of  the  envi¬ 
ronment  process.  It  will  be  instructive  to  consider  the  prob¬ 
lem  of  estimating  the  parameters  of  an  age-dependent  en¬ 
vironment  process  model.  Due  to  extensive  data  require¬ 
ments,  it  appears  that  this  extension  may  only  be  feasible 
if  the  failure  dynamics  of  the  material  under  consideration 
are  known.  Finally,  in  order  to  facilitate  a  fully-automated 
prognostic  system,  it  will  be  necessary  to  develop  more  ro¬ 
bust  techniques  for  evaluating  the  degradation  rates.  Such 
techniques  may  need  to  account  explicitly  for  the  effects  of 
external  covariates  as  noted  in  Meeker  et  al.  (2002). 
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